#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TNtuple.h"
#include "TRandom.h"
#include "TEllipse.h"
#include "TClonesArray.h"
#include "../../Objects/ObjectClass.hh"

using namespace std;

float alpha = 0.18; //Scale factor  
float PDE420nm = 0.25*alpha;
float PDE128nm = 0.1*alpha;

void read(){

  int ret = gSystem->Load( "../../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }


  int RR=20;
  int CONFIG=1;
  const int nJPMT = 19;
  float REF=0.9;

  TFile *fread=new TFile(Form("/scratch/shinichi/OutputFiles/R%d_Config%d_ref%.1f.root",RR,CONFIG,REF));
  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa[nJPMT];
  int nJPMThitsD[nJPMT];
  int nJPMThitsR[nJPMT];

  for(int iJPMT=0;iJPMT<nJPMT;iJPMT++){
    JPMTa[iJPMT]=new TClonesArray("PMT");
    tree->GetBranch(Form("JPMT%db",iJPMT))->SetAutoDelete(kFALSE);
    tree->SetBranchAddress(Form("JPMT%db",iJPMT),&JPMTa[iJPMT]);
    tree->SetBranchAddress(Form("nJPMT%dhitsD",iJPMT),&nJPMThitsD[iJPMT]);
    tree->SetBranchAddress(Form("nJPMT%dhitsR",iJPMT),&nJPMThitsR[iJPMT]);
  }

  TClonesArray *Photona=new TClonesArray("Particle");
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("Particleb",&Photona);

  int nPhoton;
  float InitXpos,InitYpos,InitZpos;
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("InitXpos",&InitXpos);
  tree->SetBranchAddress("InitYpos",&InitYpos);
  tree->SetBranchAddress("InitZpos",&InitZpos);

  TH2F *hXZInit=new TH2F("hXZInit","hXZInit",500,-250,250,640,-820,-180);
  TH2F *hXYInit=new TH2F("hXYInit","hXYInit",500,-250,250,500,-250, 250);
  TH2F *hXZDirect=new TH2F("hXZDirect","hXZDirect",50,-250,250,64,-820,-180);
  TH2F *hXYDirect=new TH2F("hXYDirect","hXYDirect",50,-250,250,50,-250, 250);
  TH2F *hXZReflect=new TH2F("hXZReflect","hXZReflect",50,-250,250,64,-820,-180);
  TH2F *hXYReflect=new TH2F("hXYReflect","hXYReflect",50,-250,250,50,-250, 250);
  TH2F *hXZHit=new TH2F("hXZHit","hXZHit",500,-250,250,640,-820,-180);
  TH2F *hXYHit=new TH2F("hXYHit","hXYHit",500,-250,250,500,-250, 250);

  TH1F *hVecX=new TH1F("hVecX","hVecX",2000,-200,200);
  TH1F *hVecY=new TH1F("hVecY","hVecY",2000,-200,200);
  TH1F *hVecZ=new TH1F("hVecZ","hVecZ",2000,-200,200);

  TH1F *hNDirect=new TH1F("hNDirect","hNDirect",32,-820,-180);
  TH1F *hNReflected=new TH1F("hNReflected","hNReflected",32,-820,-180);
  TH1F *hNAll=new TH1F("hNAll","hNAll",32,-820,-180);

  TH2F *hPMTHitDistribution=new TH2F("hPMTHitDistribution","hPMTHitDistribution",500,-250,250,500,-250, 250);
  //  TH1F *hPMTnHitDirect=new TH2F("hPMTnHitDirect","hPMTnHitDirect",);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  for(int iev=0;iev<nev;iev++){

    Photona->Clear();
    for(int iJPMT=0;iJPMT<nJPMT;iJPMT++){
      JPMTa[iJPMT]->Clear();
    }

    tree->GetEntry(iev);

    hXZInit->Fill(InitXpos,InitZpos);
    hXYInit->Fill(InitXpos,InitYpos);

    int nHitDirect    = 0;
    int nHitReflected = 0;
    int nHitAll = 0;

    for(int iJPMT=0;iJPMT<nJPMT;iJPMT++){

      nHitDirect    += nJPMThitsD[iJPMT];
      nHitReflected += nJPMThitsR[iJPMT];

      nHitAll +=  nJPMThitsD[iJPMT]+nJPMThitsR[iJPMT];

      int nhit=JPMTa[iJPMT]->GetEntriesFast();
      for(int ihit=0;ihit<nhit;ihit++){
	PMT *jpmt=(PMT*)JPMTa[iJPMT]->At(ihit);
	hPMTHitDistribution->Fill(jpmt->xpos,jpmt->ypos);
      }
    }

    hNDirect->Fill(InitZpos,nHitDirect);
    hNReflected->Fill(InitZpos,nHitReflected);
    hNAll->Fill(InitZpos,nHitAll);

    hXZDirect->Fill(InitXpos,InitZpos,nHitDirect   /(float)nPhoton);
    hXZReflect->Fill(InitXpos,InitZpos,nHitReflected/(float)nPhoton);

    int ntrack = Photona->GetEntriesFast();
    int trackid;
    float xpre,ypre,zpre;
    float xpos,ypos,zpos;

    for(int itrack=0;itrack<ntrack;itrack++){
      Particle *photon = (Particle*)Photona->At(itrack);

      trackid = photon->TrackID;
      xpre = photon->xpre;
      ypre = photon->ypre;
      zpre = photon->zpre;
      xpos = photon->xpos;
      ypos = photon->ypos;
      zpos = photon->zpos;

      if(xpre==InitXpos&&ypre==InitYpos&&zpre==InitZpos){
	float vecx = xpos-xpre; 
	float vecy = ypos-ypre; 
	float vecz = zpos-zpre; 
	hVecX->Fill( vecx );
	hVecY->Fill( vecy );
	hVecZ->Fill( vecz );
      }
    }
  
  }


  //  gStyle->SetOptStat(000000);
  TCanvas *c1=new TCanvas("c1","c1",1200,600);
  c1->Divide(2,1);
  c1->cd(1);
  hXZInit->Draw();
  c1->cd(2);
  hXYInit->Draw();


  TCanvas *c3=new TCanvas("c3","c3",1000,800);
  c3->Divide(2,1);
  c3->cd(1);
  hXZDirect->Draw("colz");
  c3->cd(2);
  hXZReflect->Draw("colz");

  TCanvas *c4=new TCanvas("c4","c4",1000,800);
  c4->Divide(2,2);
  c4->cd(1);
  hVecX->Draw();
  c4->cd(2);
  hVecY->Draw();
  c4->cd(3);
  hVecZ->Draw();

  TCanvas *c5=new TCanvas("c5","c5",800,600);
  TH1 *fr=gPad->DrawFrame(-820,0,-180,1.1);
  fr->SetXTitle("Z (mm)");
  fr->SetYTitle("Ratio");

  TLine *ll=new TLine(-820,1,-180,1);
  ll->SetLineStyle(21);
  ll->Draw();

  for(int ibin=0;ibin<hNAll->GetNbinsX();ibin++){
    int cont = hNAll->GetBinContent(ibin+1);
    if(cont==0) continue;
    hNDirect->SetBinContent(ibin+1,hNDirect->GetBinContent(ibin+1)/(float)cont);
    hNReflected->SetBinContent(ibin+1,hNReflected->GetBinContent(ibin+1)/(float)cont);
  }

  //  hNAll->Draw();
  hNReflected->SetLineColor(4);
  hNReflected->Draw("same");
  hNDirect->SetLineColor(2);
  hNDirect->Draw("same");

  TLegend *lge=new TLegend(0.65,0.7,0.85,0.83);
  lge->SetFillStyle(0);
  lge->AddEntry(hNDirect,"#color[2]{Direct}","l");
  lge->AddEntry(hNReflected,"#color[4]{Reflected}","l");
  lge->Draw();
  //  hPMTHitDistribution->Draw("colz");



}


void write(){

  int ret = gSystem->Load( "../../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  const int nEr=23;
  float Er[nEr]={1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,30.,40.,50.,70.,100.,150.,200.,300.,400.,500.,1000.}; //keV
  float Wpe = 20.0 ; //eV

  int RR=20;
  int CONFIG=1;
  const int nJPMT = 20;
  float REF=0.9;

  TFile *fread=new TFile(Form("/scratch/shinichi/OutputFiles/R%d_Config%d_ref%.1f.root",RR,CONFIG,REF));
  TTree *tree=(TTree*)fread->Get("tree");
  
  int nJPMThitsD[nJPMT];
  int nJPMThitsR[nJPMT];

  for(int iJPMT=0;iJPMT<nJPMT;iJPMT++){
    tree->SetBranchAddress(Form("nJPMT%dhitsD",iJPMT),&nJPMThitsD[iJPMT]);
    tree->SetBranchAddress(Form("nJPMT%dhitsR",iJPMT),&nJPMThitsR[iJPMT]);
  }

  int nPhoton;
  float InitXpos,InitYpos,InitZpos;
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("InitXpos",&InitXpos);
  tree->SetBranchAddress("InitYpos",&InitYpos);
  tree->SetBranchAddress("InitZpos",&InitZpos);

  TFile *fwrite=new TFile(Form("./rootfile/R%d_config%d_Alpha0.18_ref%.1f.root",RR,CONFIG,REF),"recreate");

  TNtuple *nt=new TNtuple("nt","nt","nPhoton:nHitDirect:nHitReflected:nHitAll");
  TH1F  *hNpe[nEr];
  for(int iEr=0;iEr<nEr;iEr++){
    hNpe[iEr]=new TH1F(Form("hNpe_%dkeV",(int)Er[iEr]),Form("hNpe_%dkeV",(int)Er[iEr]),5000,0,5000);
  }
  TH1F *hEffvsEr=new TH1F("hEffvsEr","hEffvsEr",1000,0,1000);
  TH2F *hRatiovsZpos=new TH2F("hRatiovsZpos","hRatiovsZpos",640,-820,-180,1000,0,1);

 
  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  for(int iev=0;iev<nev;iev++){

    tree->GetEntry(iev);

    int nHitDirect    = 0;
    int nHitReflected = 0;
    int nHitAll = 0;

    for(int iJPMT=0;iJPMT<nJPMT;iJPMT++){
      nHitDirect    += nJPMThitsD[iJPMT];
      nHitReflected += nJPMThitsR[iJPMT];

      nHitAll +=  nHitDirect;
      nHitAll +=  nHitReflected;
    }

    hRatiovsZpos->Fill(InitZpos,nHitReflected/(float)nHitAll);

    float EffDirect    = nHitDirect/(float)nPhoton;
    float EffReflected = nHitReflected/(float)nPhoton;

    for(int iEr=0;iEr<nEr;iEr++){
      float MuDirect    = Er[iEr]*1e3/Wpe * EffDirect    * PDE128nm; 
      float MuReflected = Er[iEr]*1e3/Wpe * EffReflected * PDE420nm; 

      int Npe = gRandom->Poisson(MuDirect) + gRandom->Poisson(MuReflected);
      hNpe[iEr]->Fill( Npe );
    }

    nt->Fill(nPhoton,nHitDirect,nHitReflected,nHitAll);
  }

  //---------------------
  for(int iEr=0;iEr<nEr;iEr++){ 
    float eff = hNpe[iEr]->Integral(hNpe[iEr]->FindBin(10.),hNpe[iEr]->GetNbinsX()) / (float)hNpe[iEr]->Integral();
    hEffvsEr->Fill(Er[iEr],eff);
    hEffvsEr->SetBinError( hEffvsEr->FindBin(Er[iEr]), sqrt(eff*(1-eff)/(float)hNpe[iEr]->Integral()) );
  }

  fwrite->Write();

}

Double_t function(double *x,double *par){

  //Double_t fff=-1*exp(par[0]*x[0]*x[0]+par[3]*x[0]+par[1])+par[2];
  Double_t fff1=(TMath::Erf((x[0]-par[0])/par[1])+1.0)*0.5*par[2];
  Double_t fff2=(TMath::Erf((x[0]-par[3])/par[4])+1.0)*0.5*(1.0-par[2]);

  if( (fff1+fff2)<0.0001 ) return 0 ;
  else                    return fff1+fff2;
  //  return fff;
}

void test(){

  int RR=20;
  int CONFIG = 1;
  float REF = 0.9;

  TFile *test=new TFile(Form("./rootfile/R%d_config%d_Alpha0.18_ref%.1f.root",RR,CONFIG,REF));

  //-----
  const int nEr=23;
  TH1F *hNpe[nEr];
  //  float Er[nEr]={1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,30.,40.,50.,70.,100.,200.,500.,1000.}; //keV
  float Er[nEr]={1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,30.,40.,50.,70.,100.,150.,200.,300.,400.,500.,1000.}; //keV
  for(int i=0;i<nEr;i++){
    char name[128];
    sprintf(name,"hNpe_%dkeV",(int)Er[i]);
    hNpe[i]=(TH1F*)test->Get(name);
  }

  //-----
  TH1F *heff=(TH1F*)test->Get("hEffvsEr");
  //  TF1 *func=new TF1("func","-1*exp([0]*x*x+[1])+[2]",0,1000);
  //  TF1 *func=new TF1("func","-1*exp([0]*x*x+[1]+[3]*x)+[2]",0,1000);
  TF1 *func=new TF1("func",function,0,1000,5);

  /*
  func->SetParameter(0,-0.019);
  func->SetParameter(1,0.0032);
  func->SetParameter(2,0.98);
  func->SetParameter(3,0.036);
  */
  
  func->SetParLimits(0,0,100.);
  //func->SetParameter(0,6.5);
  func->SetParLimits(1,0.,100.);
  func->SetParLimits(2,0.1,1);
  func->SetParLimits(3,10.,100.);
  func->SetParLimits(4,10.,100);


  gStyle->SetOptFit(1);

  TCanvas *c1=new TCanvas("c1","c1",800,600);
  c1->SetGrid();
  TH1 *frame1=gPad->DrawFrame(0,0,200,1);
  frame1->SetTitle(Form("Config:%d  /  #alpha = 0.18",CONFIG));
  frame1->SetXTitle("Er (keV)");
  frame1->SetYTitle("Efficiency");
  heff->SetMarkerStyle(21);
  heff->SetMarkerSize(1.);
  heff->Draw("samepe");
  heff->Fit("func","0");
  func->Draw("same");

  TCanvas *c2=new TCanvas("c2","c2",800,600);
  c2->SetGrid();
  TLegend *lg=new TLegend(0.7,0.5,0.85,0.9);
  lg->SetFillStyle(0);

  TH1 *frame2=gPad->DrawFrame(0,0,50,300);
  frame2->SetTitle(Form("Config:%d  /  #alpha = 0.18",CONFIG));
  frame2->SetXTitle("N_{pe}");
  frame2->SetYTitle("NEV");
  for(int i=0;i<9;i++){
    hNpe[2*i]->SetLineWidth(2.5);
    hNpe[2*i]->SetLineColor(6*i+51);
    hNpe[2*i]->Draw("same");
    lg->AddEntry(hNpe[2*i],Form("#color[%d]{%d keV}",6*i+51,(int)Er[2*i]),"l");
  }
  lg->Draw();


  cout<<func->Integral(0,1000)<<" "<<endl;

  // c1->SaveAs( Form("./plot20130318/R%d_config%d_Alpha0.18_ref%.1f_c1.png",RR,CONFIG,REF) );
  // c1->SaveAs( Form("./plot20130318/R%d_config%d_Alpha0.18_ref%.1f_c1.C",RR,CONFIG,REF) );
  // c2->SaveAs( Form("./plot20130318/R%d_config%d_Alpha0.18_ref%.1f_c2.png",RR,CONFIG,REF) );
  // c2->SaveAs( Form("./plot20130318/R%d_config%d_Alpha0.18_ref%.1f_c2.C",RR,CONFIG,REF) );

}
